{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sparse Gaussian Process Regression (SGPR)\n",
    "\n",
    "## Overview\n",
    "\n",
    "In this notebook, we'll overview how to use SGPR, the method of http://proceedings.mlr.press/v5/titsias09a/titsias09a.pdf in which the inducing point locations are learned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Make plots inline\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Data\n",
    "\n",
    "For this example notebook, we'll be using the `elevators` UCI dataset used in the paper. Running the next cell downloads a copy of the dataset that has already been scaled and normalized appropriately. For this notebook, we'll simply be splitting the data using the first 80% of the data as training and the last 20% as testing.\n",
    "\n",
    "**Note**: Running the next cell will attempt to download a ~400 KB dataset file to the current directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "import os.path\n",
    "from scipy.io import loadmat\n",
    "from math import floor\n",
    "\n",
    "if not os.path.isfile('elevators.mat'):\n",
    "    print('Downloading \\'elevators\\' UCI dataset...')\n",
    "    urllib.request.urlretrieve('https://drive.google.com/uc?export=download&id=1jhWL3YUHvXIaftia4qeAyDwVxo6j1alk', 'elevators.mat')\n",
    "    \n",
    "data = torch.Tensor(loadmat('elevators.mat')['data'])\n",
    "X = data[:, :-1]\n",
    "X = X - X.min(0)[0]\n",
    "X = 2 * (X / X.max(0)[0]) - 1\n",
    "y = data[:, -1]\n",
    "\n",
    "train_n = int(floor(0.8*len(X)))\n",
    "\n",
    "train_x = X[:train_n, :].contiguous().cuda()\n",
    "train_y = y[:train_n].contiguous().cuda()\n",
    "\n",
    "test_x = X[train_n:, :].contiguous().cuda()\n",
    "test_y = y[train_n:].contiguous().cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([16599, 18])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the GP Model\n",
    "\n",
    "We now define the GP model. For more details on the use of GP models, see our simpler examples. This model constructs a base scaled RBF kernel, and then simply wraps it in an `InducingPointKernel`. Other than this, everything should look the same as in the simple GP models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gpytorch.means import ConstantMean\n",
    "from gpytorch.kernels import ScaleKernel, RBFKernel, InducingPointKernel\n",
    "from gpytorch.distributions import MultivariateNormal\n",
    "\n",
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = ConstantMean()\n",
    "        self.base_covar_module = ScaleKernel(RBFKernel())\n",
    "        self.covar_module = InducingPointKernel(self.base_covar_module, inducing_points=train_x[:500, :], likelihood=likelihood)\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return MultivariateNormal(mean_x, covar_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "likelihood = gpytorch.likelihoods.GaussianLikelihood().cuda()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood).cuda()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/25 - Loss: 0.796\n",
      "Iter 2/25 - Loss: 0.786\n",
      "Iter 3/25 - Loss: 0.773\n",
      "Iter 4/25 - Loss: 0.762\n",
      "Iter 5/25 - Loss: 0.748\n",
      "Iter 6/25 - Loss: 0.735\n",
      "Iter 7/25 - Loss: 0.724\n",
      "Iter 8/25 - Loss: 0.711\n",
      "Iter 9/25 - Loss: 0.698\n",
      "Iter 10/25 - Loss: 0.685\n",
      "Iter 11/25 - Loss: 0.670\n",
      "Iter 12/25 - Loss: 0.657\n",
      "Iter 13/25 - Loss: 0.645\n",
      "Iter 14/25 - Loss: 0.631\n",
      "Iter 15/25 - Loss: 0.617\n",
      "Iter 16/25 - Loss: 0.602\n",
      "Iter 17/25 - Loss: 0.588\n",
      "Iter 18/25 - Loss: 0.574\n",
      "Iter 19/25 - Loss: 0.561\n",
      "Iter 20/25 - Loss: 0.545\n",
      "Iter 21/25 - Loss: 0.529\n",
      "Iter 22/25 - Loss: 0.515\n",
      "Iter 23/25 - Loss: 0.500\n",
      "Iter 24/25 - Loss: 0.484\n",
      "Iter 25/25 - Loss: 0.470\n",
      "CPU times: user 10.2 s, sys: 13.2 s, total: 23.4 s\n",
      "Wall time: 3.29 s\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.SGD(model.parameters(), lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "training_iterations = 25\n",
    "def train():\n",
    "    for i in range(training_iterations):\n",
    "        # Zero backprop gradients\n",
    "        optimizer.zero_grad()\n",
    "        # Get output from model\n",
    "        output = model(train_x)\n",
    "        # Calc loss and backprop derivatives\n",
    "        loss = -mll(output, train_y)\n",
    "        loss.backward()\n",
    "        print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iterations, loss.item()))\n",
    "        optimizer.step()\n",
    "        torch.cuda.empty_cache()\n",
    "        \n",
    "# See dkl_mnist.ipynb for explanation of this flag\n",
    "with gpytorch.settings.use_toeplitz(True):\n",
    "    %time train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making Predictions\n",
    "\n",
    "The next cell makes predictions with SKIP. We use the same max_root_decomposition size, and we also demonstrate increasing the max preconditioner size. Increasing the preconditioner size on this dataset is **not** necessary, but can make a big difference in final test performance, and is often preferable to increasing the number of CG iterations if you can afford the space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "likelihood.eval()\n",
    "with gpytorch.settings.max_preconditioner_size(10), torch.no_grad():\n",
    "    with gpytorch.settings.use_toeplitz(False), gpytorch.settings.max_root_decomposition_size(30), gpytorch.settings.fast_pred_var():\n",
    "        preds = model(test_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test MAE: 0.0909833088517189\n"
     ]
    }
   ],
   "source": [
    "print('Test MAE: {}'.format(torch.mean(torch.abs(preds.mean - test_y))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
